Efficient classical density-functional theories of 
rigid-molecular fluids and a simplified free energy 
functional for liquid water 



Ravishankar Sundararaman, T. A. Arias 
Department of Physics, Cornell University, Ithaca, NY 14-853, USA 



Abstract 

Classical density-functional theory provides an efficient alternative to molec- 
ular dynamics simulations for understanding the equilibrium properties of inho- 
mogeneous fluids. However, application of density- functional theory to multi- 
site molecular fluids has so far been limited by complications due to the implicit 
molecular geometry constraints on the site densities, whose resolution typically 
requires expensive Monte Carlo methods. Here, we present a general scheme 
of circumventing this so-called inversion problem: compressed representations 
of the orientation density. This approach allows us to combine the superior 
iterative convergence properties of multipole representations of the fluid con- 
figuration with the improved accuracy of site-density functionals. Next, from 
a computational perspective, we show how to extend the DFT-t— t- algebraic 
formulation of electronic density-functional theory to the classical fluid case 
and present a basis-independent discretization of our formulation for molecular 
classical density-functional theory. Finally, armed with the above general frame- 
work, we construct a simplified free-energy functional for water which captures 
the radial distributions, cavitation energies, and the linear and non-linear dielec- 
tric response of liquid water. The resulting approach will enable efficient and 
reliable first-principles studies of atomic-scale processes in contact with solution 
or other liquid environments. 



1. Introduction 

The microscopic structure of liquids plays an important role in several bio- 
logical processes and chemical systems of technological importance, and is the 
subject of continued scientific interest. Several computational techniques have 
been developed to study the bulk and inhomogeneous properties of liquids. (See 
[1] for a comprehensive review.) 

Monte Carlo calculations and molecular dynamics simulations with a sim- 
plified Hamiltonian, often composed of additive pair-potentials, are the most 
popular techniques used to compute properties of inhomogeneous liquids. How- 
ever these can be quite expensive due to the long equilibration times and exten- 
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sive phase-space sampling necessary to compute thermodynamic averages with 
sufficiently low statistical noise. 

Theories in terms of the equilibrium densities rather than individual con- 
figurations of molecules avoid this phase space sampling and hence are much 
more efhcient for the computation of equilibrium properties. Integral equation 
theories, based on approximations to the diagrammatic series of interactions, 
can be reasonably accurate but still prove relatively expensive and have only 
recently been applied to inhomogeneous systems in three dimensions [2]. 

All of the above methods require the construction of a simplified Hamilto- 
nian, usually restricted to pair potentials. Many applications, such as the deter- 
mination of chemical reaction pathways, also require estimation of free energies, 
which involves a coupling constant integration and hence incurs additional costs. 
Classical density-functional theories based on an exact variational theorem for 
the free energy of a liquid ^ avoid these restrictions, at least in principle. In 
practice, they involve directly approximating the free energy as a functional of 
the liquid density. They have the further advantage of being readily coupled to 
a quantum mechanical calculation of an electronic system within the framework 
of joint density- functional theory @], which makes quantum treatment practical 
for much larger systems than possible with ab initio molecular dynamics [5]. 

Free-energy functional approximations for fluids of spherical particles often 
employ a thermodynamic perturbation about the hard sphere fluid described 
accurately by fundamental measure theory [SI [7] . These may be extended to 
model polar fluids such as the Stockmeyer fluid jS], but the accuracy of such 
theories for real molecular fluids is not satisfactory. 

Molecular fluids are best described within the reduced interaction-site mod- 
els (RISM) [9j, which express the interactions in terms of a few sites on each 
molecule, usually on atomic centers constrained by a rigid model molecular ge- 
ometry. The free energy functional descriptions in terms of these site densities, 
however, is complicated by the molecular geometry constraints; even the ideal- 
gas free-energy is no longer expressible as an analytical closed-form functional 
of the site-densities alone. An explicit functional can be written by introduc- 
ing effective ideal-gas site potentials as auxiliary variables [TD], but this still 
requires inversion of an integral equation to obtain these potentials from the 
site densities, a problem which can be solved explicitly only in some limits such 
as reducing the molecule to a point , and requires an expensive Monte Carlo 
integration for the general case. 

The above inversion problem can effectively be avoided [12] by switching to 
the site potentials as the independent variables instead of the site densities. This 
method was applied successfully to fluids of hydrogen chloride [12] and water 
[13j in one dimension. The convergence of free energy minimization with respect 
to these independent variables turns out to be quite slow, however, particularly 
in the presence of strong electric flelds. 

This work presents a simple general scheme of choosing independent variables 
that can generate the site densities for the free-energy functional treatment of 
molecular fluids. In Section [2] we demonstrate the site-potential solution as 
a special case of this general scheme and present other representations with 



2 



better iterative convergence during free energy minimization. In section [3j we 
construct a simplified semi-empirical excess functional for liquid water which 
adequately captures the properties most critical to successful ab initio treatment 
of solvation within the framework of joint density-functional theory. Finally in 
section |4j we detail the computational implementation of the above theories 
in the open-source plane- wave density- functional theory software JDFTx [T3], 
using the basis-independent DFT-|--f algebraic formulation [TS], and present 
numerical studies of the molecular classical density-functional framework and 
the free energy functional for liquid water. 



2. Free energy of an ideal gas of rigid molecules 

The site-density-functional theory of molecular fluids is based on functional 
approximations to the in-principle exact free-energy functional 

'^[{Nair^}] = *id[{iVa}] + F,^[{N^}], (1) 

where $ is the grand free energy of the interacting fluid, $id is the exact grand 
free energy for the molecular ideal gas, Nair) are the densities of distinct sites 
(indexed by a) in the molecule, and F^x captures the effect of all the interactions 
and correlations. The equilibrium densities and free energy are obtained by 
minimizing the free energy over all allowed densities. 

The heart of the inversion problem lies in the fact that the site densities 
A^a(r) are not independent variables, but are constrained by the assumption 
of a rigid-molecular geometry. For definiteness, let the molecule geometry be 
specified by Rak , the positions of the sites for a molecule centered at the origin 
in some reference orientation. Here, a indexes distinct sites while k indexes 
multiple sites of the same type equivalent under the symmetry of the molecule 
(e.g. for a 3-site water model, a € {O, H}, k = 1 for a — O and k e {1, 2} for 
a = H.) 

2.1. Treatment o J site- density constraints 

The inversion problem in the original approach of (TUl E] , which includes 
ideal-gas effective potentials ipai^ as auxiliary independent variables in addi- 
tion to the site densities, is avoided in [12] by switching to ipair) as the sole 
independent variables. The site densities and ideal-gas free energy in the pres- 
ence of external site potentials Va and chemical potentials fj,a are then expressed 
in terms of the ipa (t) using 



<fid[{Mr)}] = f^^^^Wa}] + J2 / driV„(f)(K.(rO - - ^,(f)) (2) 
^ -TV.efT J l[e-^-^^^-^^/^dr^^ks{{r^k}) (3) 

a.k 
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Here, the reference density A^ref sets the zero of chemical potential and the con- 
straint function s{{fak}) picks out configurations {^0,^} which satisfy the rigid 
molecule geometry (i.e. equivalent to {Rak} under rotations and translations). 

Employing a spherical harmonic expansion of the constraint function, |12) 
and [T3] specialize ([s]) for diatomic and triatomic molecules respectively. How- 
ever, that expansion also becomes computationally challenging as one moves to 
calculations without planar symmetry. Instead, we transform ([3| to 

= -iV,efT / ^ n e-^"(^"+"°«-)/^ (5) 

where uj e 5*0(3) is a rotation and a; o _R is the result of rotating vector R 
by w, and we directly discretize the integral over orientations as described in 
[Appendix A for practical calculations in three dimensions. 



It is instructive to further transform the above equations to 

= ^P-(^ (log ^ - 1) + E / dr7Va(rl(F„(rl - ^^^) (6) 
(^1 = E / ^P'^ (^^ - ^ ° ^-k) (7) 



with 



P.(rO=7V„fne-^°^'^+"°'^°'=^/^- (8) 



a.k 



Here, Pui{r) represents the probability of finding a molecule centered at location 
r with orientation lu. For an ideal molecular gas, Puj{r) is simply a product of 
Boltzmann factors for each site given that V'a (r) ideal-gas effective potentials 
since they equal Va (r) — when $id is minimized. 

Note that, given the explicit expressions for the ideal-gas free energy ^ and 
site densities ([t]) , (^ is a natural choice for the independent variables for un- 



constrained free-energy minimization. Section 4.2 demonstrates that conjugate 
gradients minimization over p^^ (r) as the independent variables converges much 
faster than minimization over the {ipair}}- A potential disadvantage of using 
PojiT) is the increased memory requirement, but practical calculations of reason- 
able size are possible using the efficient orientation quadratures of |Appendix A[ 
Moreover the superior convergence properties can be retained, while mitigating 
the memory requirements, by switching to compressed multipole representations 
oipui{r), as we now discuss. 

2.2. Representations of the Orientation Density 

Our first task in this development is to demonstrate that minimizing free 
energy functionals over p^^{r) yields the same results as minimizing over {ipaif)}- 
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To demonstrate this we employ a constrained search procedure, to find 
$ = min ($id[p.(r)] + F,^[{N^}]) 

= mm (t f ^p.(r01og^+FM-e.[{A^4]), 

which follows because all the terms in $id are explicit site-density function- 
al except for the molecular ideal gas entropy (S'id) contribution separated out 
above. Next, the minimization over all Puj{r) can be performed by minimizing 
over those PLj{r) that yield a specific set of site densities {Naif)}, and then 
minimizing over all {Nair)}, 

min ( min T / ^p^(r) log ^ + Fid.ex[{^a}]V 

Finally, the inner, constrained minimization over Pu}{r) that lead to given site 
densities can be performed explicitly by introducing Lagrange multipliers Tpa {"r) 
for each Na{f) constraint. It is straightforward to verify that the Euler-Lagrange 
equation for that extremization is precisely ([s]) , so that the result of free energy 
minimization over (r) is exactly the same as the ideal-gas effective potential 
methods of [ini[I2]- 

To generalize this approach, we note that the exact equivalence between 
minimization over Pu>{r) and minimization over {i^ai/)} holds only when the 
external potential takes the form of external site potentials Va (r) . In principle, 
we could go beyond the reduced-interaction site model and consider arbitrary 
orientation dependent external potentials (r) (of which site potentials Va (f) 
are a special case). From this perspective, the minimization over {V'a(r)} '^^^ be 
reinterpreted as a minimization over only those p^j (f) that maximize the molec- 
ular ideal gas entropy Sid[Pui{'i^] subject to site-density constraints {Na{r)} 
(whose Lagrange-multiplier constraints become the site potentials.) The vari- 
ational principle implies that this procedure will always result in a free-energy 
greater than or equal to direct, unconstrained minimization over Pui{'r), with 
equality guaranteed only when the external orientation potential VL(r) can be 
reduced to site potentials Vair)- 

These considerations lead to the perspective of the {V'a(^)} a-s a compressed 
representation of Puj{r), with decompression carried out by maximizing the 
entropy subject to constraints for which the {V'a(^0} are Lagrange multipli- 
ers. From the most general perspective, then, any set of functional constraints 
{Xi = Xi[pi^{f)]} corresponds to a maximum-entropy compressed representa- 
tion of Puj{r), where the independent variables Xi for the free-energy functional 
minimization are the Lagrange multipliers for the corresponding Xi constraint 
in the maximization of Sid[Puii^]- Specifically, 

$ = min (c&id K(rO[x.]] + ^^ex [{iV„K(r)[x.]]}] , ) (9) 

iXi} 



5 



where {r) [xi] is the solution of 

(^r5id[p.(r-)] + Y.^MPu.m - X^)x}j = 0. (10) 

Here, $id[Pw(^)] and Na[pu{r)\ are given by ([6| and ([t]) respectively. Note that 
i typically includes a continuous index such as r, and then denotes the 
corresponding integrals. 

From this new perspective, picking Xi[p^{r)\ ~ Na[Pi^{r)\ yields the ideal- 
gas site-potential representation with Xi = '/'a(^ as the independent variables 
and Puii'r) given by Similarly, picking Xi[pi^{r)\ = Puii'r) yields the trivial 
self-representation, with Pujij) as the independent variables. As shown earlier, 
both these representations are exact when the external potentials are site poten- 
tials, while the former is a variational approximation to the latter in the most 
general case of orientation potentials. 

The advantage of this general framework is that we can develop new, phys- 
ically motivated representations which then are guaranteed to be variational 
approximations. Of particular interest are representations based on multipole 
probability densities 

ML,,n. (r-) [p^ (f)] = I du;p^ (r)i^4^,„, {uj) , (11) 

where D^^^^{uj) are the Wigner D-matrices [H] (irreducible matrix represen- 
tations of SO{3)). The Lagrange-multiplier independent variables /J.4,im2(r) 
resulting from this choice then generate the orientation probability 

pUrl[t^L.mArl]=Nr.iYl ft expf- ^-^"-^"^ii^-^--^"^ ). (12) 



J mi,m2=-J 



By the completeness of the Dif^^^^ on 50(3), this representation is exact if all 
components j — > c» are included. In practice, we truncate the expansion at 
finite 

We find below that including terms up to j = 1 is sufficient for many prac- 
tical problems, particularly when the external potential is dominated by strong 
electric fields. We choose to label the corresponding independent variables for 
this truncation as fi{r) for j — and e(r) for j — 1, because they correspond 
to the ideal-gas effective local chemical potential and local electric field (up to 



factors of T and the molecule's dipole moment). Section 4.2 below compares 
the accuracy and convergence properties of this {/i, e| representation to those 
of the site-potential {{ipa}) representation and the self-representation (pu)- 



^ This expansion in j is different from the spherical harmonic expansion of for tri- 

atomic molecules introduced in I13j . In particular, truncating expansion | |12| l at j = 1 retains 
the exact nonlinear dielectric response for axisymmetric molecules, whereas the corresponding 
truncation in |13| would incur a 20% error in the O(E^) term of t(E) at ambient conditions. 
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Finally, we would like to point out that this general perspective opens up a 
promising avenue for excess functional development. Our framework enables the 
computation of site densities and multipole densities irrespective of the inde- 
pendent variables used for minimization, which facilitates the generalization of 
site-density excess functionals Fex[{-^a}] to combined site- multipole functionals 
-Fgx[{^q}, {A^fijiima}] eyeu to full orientation density functionals Fex[Puj]- In 
particular, it should now be possible to combine the best features of site-density 
functionals, which better capture short-ranged correlations, with those of multi- 
pole functionals, which allow for analytically derivable long-range correlations. 



3. Excess functionals 



So far we have focused on accurate and efficient representations of the ideal 
gas of rigid molecules. These need to be combined with good approximations 
for the excess functional Fex[{-^a}] to obtain a practicable theory for inhomo- 
geneous liquids. 

3.1. Excess functionals for model fluids 

The fluid of hard spheres has been studied extensively theoretically as well 
as with computer simulations. Within classical density-functional theory, it is 
accurately described by Rosenfeld's fundamental measure theory [5], which sat- 
isfies several rigorous conditions such as reducing to the exact Percus functional 
in the inhomogeneous one dimensional limit |17j and reproducing the Percus- 
Yevick pair correlations |18j in the bulk three dimensional limit. 

There are several variants of the fundamental measure theory functional 
corresponding to different bulk equations of state and regularizations for the 
zero-dimensional limit. (See |7] for a detailed review.) The excess functional i^ex 
for the highly accurate 'White Bear mark IF variant [11] based on the Carnahan- 
Starling equation of state for the bulk hard sphere fluid [20|, including tensor 
regularizations due to Tarazona |21| . is 



'^m[N]=T / dr 



/ ^lolog:; h/2(n3) h 

1-^3 1 - "3 



fain^)- 



3n2|nt,2p -I- 9 (jiy 



2 • ■ nv2 



2 



with 

/2(n3)-l + 
/3("3) = 1- 



247r(l - 713)2 

713(2 - ns) + 2(1 - Us) log(l - 713) 



(13) 



and 



2n3 - 3?7§ + 2nl + 2(1 - 773)^ log(l - 773) 



3nl 



where the n^'s are scalar (7 = 0,1,2,3), vector (7 = 7;1, i;2) and rank-2 tensor 
(z = m2) weighted densities defined as ni{'r) = Wi * N = J dr'wi{'r — r')N{f') 
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for hard sphere density N{f). The weight functions Wi are spherical measures 
of various dimensions (volume, surface etc.) given by 



wo{r) = ^(-Rhs 



r)/(47rr^) 
r)/(47rr) 
W2(r) = 5(i?HS - r) 
wsir) = 0{Riis - r) 



m„i(f) 



-5(i?HS - r) 



rr 1 =!\ 
/ 



S{Rm " r) 



(14) 



The hard sphere fluid also serves as an excellent reference for perturbation 
theory for other model systems. For example, the pair-potential for the Lennard- 
Jones fluid 



f/Lj = 4e 



(- 



(15) 



with energy scale parameter e and range parameter a is often split into repulsive 
and attractive parts [H] as 



UA{r) = 



e + 4e 
0, 



4e 



(7)" -(f)' 



, r < 

r > 7}l^u 

r < 
r > 



(16) 
(17) 



The free energy functional for this fluid can be approximated by treating the 
fluid interacting with Uji{r) alone using fundamental measure theory, typically 
with a hard sphere radius i?HS = f/2, and then accounting for the effects of 
UA{r) perturbatively. Mean field perturbation then leads to the excess func- 
tional 

Fi^^^iNir}] « $Hs[7V] + \ Uf I dr'7V(r)C/^(|f - r'|)7V(r'), (18) 



and several beyond-mean-field approaches have been developed to improve upon 
this starting point. 

Of particular interest is the recent approach of Peng and Yu [2^ to recast 
the mean-field term into a nonlinear weighted-density form 



j,(MWF)[^(^-.)] « <i>Hs[A^] + / dTN{T)All,{wA * N) 



(19) 
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with the mean-field weight function set to the normafized perturbation potential 



Here, ^att(^) = ^lj(^) - ^hs(^) is the difference between the Helmholtz 
energy per particle for the uniform Lennard- Jones fluid and the uniform hard 
sphere fluid at the same bulk density N. Peng and Yu demonstrate that this 
functional does an excellent job of reproducing the inhomogeneous density pro- 
files and vapor-liquid interface energies in comparison to Monte Carlo simula- 
tions of the Lennard-Jones fluid. 

3.2. Excess functional for liquid water 

The situation for a polar molecular fluid such as water is much more com- 
plicated than the model fluids mentioned above. Most approaches to the excess 
free energy of inhomogeneous water |M1 HTJ |T3] are constructed to reproduce 
the pair-correlations in the uniform fluid limit obtained by computer simulations 
or from neutron-scattering data. They can be reasonably accurate for modest 
inhomogeneities, but their practicality is limited as they are tied to the temper- 
ature and pressure of the simulation/experiment data that they are based on, 
and usually lack a simple analytic formulation. 

An alternate strategy is based on identifying a simple model Hamiltonian 
for which an approximate analytic free energy functional is readily formulated, 
and then constraining the parameters of the model Hamiltonian to the bulk 
properties of the fluid, such as the equation of state. Wertheim's thermody- 
namic perturbation theory [25 is a useful framework for generating free energy 
functionals; one class of Hamiltonians considered for water within this frame- 
work is based on tetrahedral association sites for hydrogen bonds HHj , but these 
models are yet to successfully predict the quantities relevant to solvation such 
as pair correlations, cavitation energies and dielectric response, partly due to 
the relative complexity of the model Hamiltonian. 

Recently, we proposed an alternate model Hamiltonian 27] based on cap- 
turing the effects of the empty space in the tetrahedral hydrogen bond network 
by attaching 'void' spheres to the molecule in the directions conjugate to the 
tetrahedral hydrogen-bond directions. The bonding constraints in the result- 
ing rigid trimers of hard spheres was also treated using Wertheim perturbation 
theory, but the relative simplicity of that model enabled an accurate free en- 
ergy functional description of the inhomogeneous fluid capable of predicting the 
aforementioned quantities relevant for solvation. 

This 'bonded- voids' free energy functional for water is adequately accurate 
for cavitation energies, dielectric response and the height and particle content 
of the first peak in the pair correlation. However, the secondary peaks in its 
pair correlation occur at the characteristic distances for a close-packed hard 
sphere fluid rather than for a tetrahedrally-bonded one. Evidently the cavitation 
energies are not sensitive to this deficiency in the secondary structure of the pair 




(20) 
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correlation; the height of the first peak and the exclusion volume (location of 
pole) in the equation of state are the important factors, which are captured 
correctly by the bonded void spheres ansatz. 

Here, we present a simplified free energy functional for water which retains 
only the critical features of the bonded- voids model ^7\, while eliminating the 
complexity of Wertheim perturbation. This functional employs a hard sphere 
reference with a weighted density term constrained to reproduce the equation of 
state in the spirit of the approach of l23] for the Lennard- Jones fluid. Due to the 
polar nature, we need to distinguish between short-ranged orientation-averaged 
interactions with a tail similar to the Lennard- Jones pair potential and 
long-range orientation-dependent interactions with ar~^ tail between individual 
charged sites resulting in r^^ for neutral molecules with a net dipole moment. 

We deal with the long range orientation-dependent part by taking advantage 
of the rigid molecule site- model capability developed in section |2] In particular, 
we adopt the molecule geometry and site charges of the popular SPC/E pair 
potential model [28) for molecular dynamics simulations of water, which consists 
of an O site with charge Zq = -1-0.8476 and two H sites with charge Zh = 
—0.4238 e" in a bent geometry with an 0-H distance of 1 A and a tetrahedral 
H-O-H angle (cos-i(-l/3) ~ 109.5°). 

For the shorter-ranged orientation-dependent part, we assume a Lennard- 
Jones interaction between the O-sites since it has the correct r^^ tail. We 
arrive at the excess functional ansatz 

^ex^''[^o(r),iVH(r)] « $Hs[A^o] + J drNo{r^A^tf{wA * No) 
, MT) ^ 



, 2 ^ df / dr'7V„(r)i^(|r-r'|)7V^(r'), (21) 

a.Pe{0,H} 

by adding a long-range polar correction (third term) to the Lennard- Jones func- 
tional of [23] (first two terms) . The following paragraphs specify the Helmholtz 
energy function A^^l'^{N), the dipole correlation factor A^{T) and the modified 



Coulomb kernel K{r). We shall refer to this excess functional (21) as 'scalar- 
EOS' because the excess free energy density is attributed to the scalar moment 
of the orientation density and is constrained to the equation of state. 

<i>HS is the White Bear mark II fundamental theory functional. 



13 1, for a fluid of hard spheres of radius i?HS- The second weighted 



In ([21 
given by ( 

density term employs the mean-field weight function wa (r) given by ( 20 ) with 
<J = 2Rm- 

The third term of pTj ) is the mean-field electrostatic interaction between the 
charge-site densities scaled by a dipole-correlation factor A^{T). Following [13], 
the Coulomb kernel K{r) is cutoff at high frequencies as 



47r 



(22) 

with Gc = 0.33 bohr~^, and the dipole correlation factor is chosen to reproduce 
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the bulk linear dielectric constant. Without the correlation factor, i.e. with 
Af^ = 1, the SPC/E geometry would yield a dielectric constant of 19.7 at ambient 
conditions instead of the experimental value of 78.4. The single parameter fit 

reproduces the bulk linear dielectric constant over the entire liquid phase with 
a relative RMS error ^1%. 

Next, we constrain F^^'-' to reproduce the correct Helmholtz energy den- 
sity for the uniform fluid of molecular density iV, which may be obtained by 



integrating the equation of state (ji{N,T)). Note that the third term of (21 1 
does not contribute to the uniform fluid free energy, and hence ^^t*^ must be 
the difference between the per-molecule Helmholtz free energy in water and the 
hard sphere fluid. Using the Jefferey- Austin equation of state [2^] for water, 
this constrains 

A2f(iV)^^log^3^-Kw + ^*T)iV 

_ 2Tr(T) '-±^ log »" + "HBe-HB/T 

^ ^'^ + Ciexp(^^^ ^ ^o + n^B 

up to a temperature-dependent constant which is absorbed into the arbitrary 



reference for the chemical potential /i. The first two lines of (24) represent 
the free energy density corresponding to the excess pressure for liquid water 
as parametrized in |29j by fits to experimental data for bulk liquid water, and 
the definitions of the numerous constants and functions of temperature may be 
found thereinj^ The last line of ( 24 1 subtracts the uniform fluid per-particle free 
energy corresponding to $hs given by (13), with Vhs = 47ri?^g/3. 



Now, \2\\ is completely specified except for the value of the hard sphere 
radius -Rhs- Unlike the Lennard- Jones fiuid, there is no prescribed pair potential 
from which it may be derived. We require that calculations with the excess 



functional (21) result in the surface-energy of the planar water liquid- vapor 
interface in agreement with the experimental surface tension of 72.0 x 10^"^ N/m 
at ambient temperature 298 K, and obtain 

i?HS = 1-36 A (25) 



The details of the planar interface calculation are presented in Section [4~1] and 
tests of the accuracy of the scalar-EOS functional for inhomogeneous liquid 
water are in Section l473l 



^ Note that the constants hsted in 1291 are in SI/CGS units, and should be converted to 
atomic units (with kg = 1) before substitution in 12411. 



11 



4. Results 

The efRcient rigid-molecular ideal gas representations of section [2] combined 
with the excess functional for water from section 3.2 forms a practical theory 
of inhomogeneous liquid water as we show below. We use this system to study 
the convergence properties of the various molecular ideal gas representations 
in section |4.2[ and then test the accuracy of the scalar-EOS water functional 
against experiment and molecular dynamics simulations in section [4. 3| 



Discretization 

The free energy functional approximations presented here involve integrals 
over space and orientations, which must all be discretized in a practical calcu- 
lation. The discretization of three dimensional space may be performed in a 
variety of bases including plane-waves, wavelets and specialized bases such as 
planar and radial one dimensional grids for high symmetry cases. 

We present the details of the numerical formulation of the free energy func- 
tional for rigid-molecular liquids using the basis-independent algebraic formu- 
lation developed for electronic density-functional theory [T5]. Within this for- 
mulation, the physics is expressed in terms of a handful of abstract operators 
independent of the basis, while the implementation of these operators in code 
is basis dependent. This allows for the same top-level physics code to be used 
with multiple basis sets with no modification. A three-dimensional plane-wave 
basis implementation of the fluid framework and excess functionals (using the 
notation and operators described below) is distributed with the open-source 
electronic density- functional theory software JDFTx [14], which specializes in 
solvated ab initio calculations. An analogous code base for high-symmetry one- 
dimensional basis sets, suitable for development and testing of new fluid func- 
tionals, is distributed as a sub-project of JDFTx 

Here, we briefly introduce the notation and operators required for classical 
density-functional theory; see [15] for a detailed description. A function of space 
/(r) is expanded in terms of basis functions {bi{r)} with coefficients /; (often 
written as a vector /) i.e. f{f) = fibi{r). 

The overlap of two functions /{r) and g{r) is 




which defines the basis overlap matrix O (which would be diagonal for orthogo- 
nal basis sets) . Similarly, any linear operator reduces to a matrix. For example, 
/ dff*{r)V'^g{r) — P Cg defines the Laplacian matrix dj — J drb*{r)V'^bj{r). 

The density functionals also involve integrals over nonlinear functions which 
of course cannot be reduced to basis-space matrices like the linear operators con- 
sidered above. Consequently, the basis sets are accompanied by a quadrature 
grid consisting of a set of nodes {rj} over which integration of nonlinear func- 
tions is performed. A function f{r) sampled on this quadrature grid fj — f{rj) 
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is denoted simply by the vector /. This introduces the Hnear basis-to-real space 
operator I defined by / = 1/ with matrix elements Iji = bi{fj), and the 
real-to-basis space operator, J — 2^icflj^ Armed with these operators, we can 
discretize the commonly encountered integral / drf{'r)A{g{r)) = pOJA{Ig) = 
f^J'^OJ'A{g) where A is some nonlinear function (which operates element-wise 
on vectors). 

In the particular case of plane- wave basis on a periodic unit cell, the quadra- 
ture grid fj is a uniform parallelepiped mesh, the basis functions are e"*'^ '" for 
reciprocal lattice vectors G, and the operators I and are implemented as 
Fast Fourier Transforms. O is the scalar matrix il, and C is the diagonal ma- 
trix — ri|Gp, where is the unit cell volume. For a detailed specification of 
these operators, see [15] for the three-dimensional plane-wave basis, [3T] for a 



multi-resolution (wavelet) basis, and Appendix B for the planar, cylindrical 
and spherical one-dimensional grids. 

In fact, the six operators introduced above (counting hermitian adjoints 
separately) are the only ones required for electronic density functional theory in 
the local density approximation (LDA) . The advantage of writing code in this 
framework is that implementing a new basis only requires reimplementing the 
small number of these operators. 

To express the classical density functionals, we need to introduce two addi- 
tional operators. Firstly, the computation of weighted densities involves con- 
volutions /i(r) = J d'r'f{f— 'r')g{r'), which may be discretized using a basis 
dependent tensor Cfj to hk = J2i j^ijfidj^ which we also denote hy h = f * g 
for brevity. Integrating the defining relation multiplied by basis functions, we 
see that the convolution tensor elements must be 

= J2(0-%i J drj dr^bnmif- r')6,(r-). (27) 

is symmetric under i j when the space is translationally invariant, and 
reduces to the element-wise multiply C^j ~ flSkiSkj for the plane-wave basis, as 
is well known. 

Secondly, the rigid molecule formalism of section [2] requires sampling func- 
tions with arbitrary displacements in order to generate orientation densities 
from the effective site potentials, and to generate the site densities from the ori- 
entation densities. This requires the inclusion of a translation operator defined 
by Tsf{r) = f{f+d) to our toolkit. This may be discretized to gi — '^j{Ts)ij fj 

where / and g are the discretizations of /(r) and f{r + a) respectively. The 
natural translation operator for a given basis set obtained by integrating the 



J = X^^ is the natural choice when the number of basis functions equals the number 
of quadrature grid points, which is the case for the plane-wave basis for example. When the 
number of grid points exceeds the number of basis functions, one possibility is to use the 
left-inverse as indicated so that = 1, although this is not necessary. 
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definition multiplied by basis functions is 



(Ta)., = J2(^''y"^ J <ir^mb,ir + a), (28) 

k 

and satisfies 7^ = T-s by definition. In the plane-wave basis, this operator 

takes the diagonal form {Ts)ij — 6ije~'-^'''^ . 

However, this 'Fourier' translation operator introduces severe ringing in func- 
tions that have components that extend up to the Nyquist frequency. This 
can be quite problematic for the classical density-functional theory of rigid 
molecules, particularly since positive functions can ring negative on translation, 
leading to regions of negative site densities even when p,^ > 0. The free energy 
functionals evaluated for negative site densities can be unphysically favorable 
which encourages further ringing, resulting in a numerical divergence]^ 

We remedy this by using inexact translation operators which have the prop- 
erty that they map the set of functions with all-non-negative samples on the 
quadrature grid onto itself. The action of the translation operator on the quadra- 
ture grid Ss = TTsSf can be viewed as sampling the function on the grid with 
displacement a. The natural translation operator for the plane-wave basis cor- 
responds to a sampling operator S based on Fourier interpolation. Amongst 
the piece-wise polynomial spline interpolations, only the constant spline (pick 
nearest neighbor) and linear spline (linear interpolation in each cell) guarantee 
non-negative results for a non-negative sample set; we denote the corresponding 
approximate sampling operators by S'~^ and respectively. 

The discretization of spatial integrals in the rigid-molecule classical density 
functional framework can be achieved using the above operators; the final in- 
gredient is the discretization of the orientation integrals. We achieve this using 
a quadrature rule directly on SO{'i)/G, where G is the symmetry group of the 
fluid molecule, so that 

^f{u^)^Y.WJ{u,) (29) 



with a finite set of orientations uji and weights Wi. Appendix A| describes var 



ious methods for the generation of quadrature rules on 5'0(3)/Z„ ranging from 
outer product quadratures on Euler angles to uniform sampling sets based on 
platonic solid rotation groups. Section [4. 2| explores the convergence of the orien- 
tation integrals with quadrature for the scalar-EOS water functional (symmetry 



group Z2), and the list of explored quadratures is summarized in Table A.l 

We can now discretize the molecular ideal gas free energy ^ given the 
orientation density pi^. on the quadrature grid for each discrete orientation and 



* In principle, the contributions to the free energy from regions of negative site densities 
could be redefined to zero. However, this results in a highly non-analytic energy landscape 
with extremely poor convergence for minimization algorithms 
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the site densities Na in the chosen basis set, as 

$id = tVOJY^W.p^^ {log j^-l]+Y. ^a'^(^" - ^"1) (30) 

Note that all unary real functions are understood to operate element-wise on 
vectors on the quadrature grid, unless otherwise specified. 

The expression of the orientation density on the quadrature grid in terms 
of the independent variables for minimization depends on the chosen represen- 
tation. In the self representation, the independent variables are p^^ in basis 
space and therefore p^o^ = Ipui ■ The independent variables in the site-potential 
representation are i/jq, and the orientation density is generated using as 



\ a,k 

= ^icf exp ^ S^^^j^^I4>^ , (31) 

where the latter expression with an approximate sampling operator S is used 
in practice. In the multipole representation, the independent variables are 
Mmim2('^ for |r7ii|,|m2| < j < imax and the orientation density is generated 
using (12) as 

w +J- \ 

= ^refexp I — -Clmim2('^«)/^mim2 L (^2) 
y j=0 mi,m2=-j j 

which simplifies for jmax = 1 in terms of independent variables /2 and e to 

-X f/i + (Wi o z) • ? j 
P^^ = iV,ef exp ^ . (33) 

Finally, the site densities are generated from the orientation density by a 
discretization of (lil given by iV„ = 4^r2("*\ with = -TX'^OJY.i ^iVu.^ 
and p^^ given by pTj), so thai[^ 

7V„ = Diag(^tOi)-i Vl^, ^■mg{J^O\)p^^. (34) 

i k 

For the translationally invariant plane- wave basis set, the above expression is 
equivalent to — X^i X^fc '^-cj ofl ^P'^i' ^^'^ intuitive discretization of (FtI), 



SThis is derived from = / ^dtp = ltC'J'Diag(|^)d</), which leads to Diag(J'tC)i) |n = 
Here, Diag(2:) is the diagonal operator with the elements of x on its diagonal, i.e. 
Diag(z)?/ = Diag(j;)x is the element-wise multiplication of x and y. 
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and this holds approximately for other three-dimensional basis sets. However, 



( 34 ) holds even when Ss is generalized to a non- uniform translation Ss{p) , which 
is required for the reduced-dimensionality basis sets of [Appendix B[ 

Movin g on to excess functionals, the hard sphere excess free energy $hs[-^] 



given by (13 1 is discretized by replacing / df — )■ VOJ and computing the 
integrand element-wise on the quadrature grid, where the weighted densities are 
computed from convolutions in basis space Ui = T{'Wi * N). These convolutions 
may be computed efficiently in the plane-wave basis by multiplying with the 



analytic Fourier transforms of the weight functions (14), but in other bases, 
they should be computed with specialized routines that take advantage of the 
finite range of the weight functions. (See [7] for examples.) The excess free 



energy of the Lennard- Jones fluid (53], given by (|19|), discretizes to Fi^^^' 



+ N^OJAl^^{I{wA * N)). The convolution WA* is trivial in the plane- 
wave basis, but may require specialized routines in other basis sets due to the 
polynomial tail of the Lennard- Jones weight functionj^ 

Finally, the scalar-EOS excess functional for water ( plj ) is discretized to 

F,^^o ^ $Hs[iVo] + iV^O^^f (I(7i^ * No)) 

+ ^4^ Z^Zp{wK*Nc.)^d{-A^C~^)d{wK*Np). (35) 



Here, the high-frequency cutoff Coulomb Kernel (22) has been rewritten in 
terms of the bare Coulomb kernel (— 47r£^^)(5 computed by solving the Poisson 
equatiorj^ by introducing the site-charge kernel 



zI;K(G) = l/yi+ — j . (36) 

This modification has no effect for the plane-wave basis, but is important for 
other basis sets because it decomposes the long-range convolution with K{G) 



into a short-ranged convolution ((36) is confined exponentially in real space), 
and the solution of Poisson equation which is a standard operation in any basis 
set [151131]. 



4-.2. Convergence 

Section |4.1| presented the discretization of the general rigid-molecular ideal 
gas framework of section [2] with various choices for the independent variables. 



° For example, in wavelet bases, this may be performed by decomposition into a finite- 
ranged part treated at all grid levels, and a bandwidth-limited long-range part performed 
using the Fourier method on the coarsest grid. 

O IS the overlap operator with the null-space of C projected out, and is understood 
to be the inverse of C in orthogonal complement of the null-space with zero projection in the 
null-space. See |15| for details. 
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and excess functionals including the scalar-EOS functional for liquid water con- 



structed in section 3.2 Next, we briefly discuss the minimization of the liquid 
free energy given a set of external potentials, compare the performance of the 
different choices of independent variables, and explore the accuracy of the dis- 
cretization of the orientation integrals. 

The free energy of the fluid for a particular excess functional and choice of 
independent variables is expressed in the basis-independent algebraic formula- 



tion of [15], including the operators introduced in section 4.1 The gradient 



of the free energy with respect to the independent variables may be derived 
in the same notation in a straightforward manner as shown in |15| . and the 
computational cost for evaluating the gradient is comparable to that for the 
free energy. We can therefore minimize the free energy functional to flnd the 
equilibrium configuration of the fluid directly using the non-linear conjugate 
gradients method |32j. 

First, we compare the convergence of the conjugate gradients method for 
different choices of independent variables. For the remainder of this section, we 
work with the scalar-EOS functional for water at a temperature of 298 K in 
the three-dimensional plane-wave basis set, and perform all calculations using 
JDFTx [M] . We focus on two physical systems which capture different extremes 
of the typical external potentials encountered in ab-initio solvation: water sur- 
rounding a hard sphere, and water in a parallel plate capacitor with a strong 
electric field (in the dielectric saturation limit). 

The hard sphere system consists of an external potential Voir) — Vo9{R — 
\r\) which excludes the O sites of water from a sphere of radius i?, with no 
potential on the H sites (Vff(r^ = 0). We pick i? = 4 A, a reasonable size 
for the region excluded by a molecule solvated in water, and Vq = 1 (« 
27.2 eV) which is sufficient to completely exclude the liquid from that region. 
The calculations are performed in a cubic unit cell of side 32 bohrs (« 17 A) 
with a 128 x 128 x 128 fast Fourier transform (FFT) grid; the grid spacing of 
0.25 bohrs corresponds roughly to the charge density grid of a typical electronic 
density-functional theory calculation at a wave-function kinetic energy cutoff of 
20 Eh. 

The parallel plate capacitor system consists of two plates 112 bohrs apart, 
with an external potential corresponding to an applied electric field of Eq = 
1 V/A (10^° V/m), which is typical in the first solvation shell of a polar molecule, 
and corresponds to a regime of strongly non-linear dielectric response. (See Fig- 
ure [sj) Repulsive potentials of strength 1 Eh on both the O and H sites confine 
the fluid to the region between the capacitor plates. The calculation is performed 
in a periodic cell of length 256 bohrs containing two capacitors back-to-back so 
that the cell has no net dipole, and is sampled using a one-dimensional FFT 
grid with 4096 points. The transverse dimensions are translationally invariant, 
and the free energies reported are per bohr^ transverse area. 

Figure [T] shows the convergence of the Polak-Ribiere variant of the nonlin- 
ear conjugate gradients algorithm [33J for the hard-sphere and capacitor systems 
with the site-potential {{tpa}), j — ^ truncated multipole ({/x, e}) and self {{puj}) 
representations of the orientation density as independent variables. The initial 
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Iterations Iterations 

(a) 4 A hard sphere (b) Capacitor with Eq — I V/A 

Figure 1: Convergence of conjugate gradients minimization of the free energy of scalar-EOS 
water (a) surrounding a 4 A hard sphere, and (b) in a paraUel plate capacitor with exter- 
nally applied field strength Eq = 1 V/A, for different independent variables. The difference 
of free energy from the final equilibrium value as a function of iteration count is shown on a 
logarithmic scale for the self representation {ptj} (solid red line), the site-potential represen- 
tation {tpa} (blue dashed line), and the multipole representation {^, e} (thicker green dotted 
line). The fainter green dotted line is the difference of the free energy in the {/J, ej- repre- 
sentation from the converged value within that representation (which is variationally higher 
than the equilibrium value). Note the rapid exponential convergence in the self and multipole 
representations, compared to the site potential representation. 
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(a) 4 A hard sphere (b) Capacitor with Eq = 1 V/A 

Figure 2: Convergence of free energy with orientation quadrature for the two systems con- 
sidered in Figure [T] The orientation quadratures studied are listed in Table |A.1[ and the 
free energy at the Euler(12) quadrature (which has 3456 nodes on SO{3) /Z2) is used as the 
reference in computing relative errors for all the smaller quadratures. Note that the error due 
to the orientation quadrature plateaus at jmax ~ 7 for the sphere geometry, and at jmax ~ 10 
for the strong-field capacitor; these would therefore be reasonable choices in ab initio solvation 
calculations for non-polar and strongly-polar molecules respectively. 



guess in each case corresponds to a uniform bulk density of water in the allowed 
regions and no density in the disallowed regions, with a uniform orientation 
distribution for the sphere geometry, and a dipolar orientation distribution cor- 
responding to bulk linear dielectric response for the capacitor geometry. The 
7-design quadrature with 96 nodes on SO{3)/Z2 (see Table A.l) was used for 
orientation sampling. 

The self representation ({p^}) exhibits the best exponential convergence, 
and is the method of choice when it is practical to store the orientation density. 
The multipole representation ({^, e}) also converges quite rapidly, but it is a 
variational approximation and will result in a higher free energy than that in 
{pui}- Note that for a typical molecule cavity formation (the hard sphere case), 
the error in the {fi, e}-representation is ^ 4 x 10^^ Eh or ^ 0.03 kcal/mol, which 
is negligible in the computation of solvation energies. Likewise, the relative er- 
ror in the free energy of the strong-field capacitor corresponds to an error of 
less than ^ 1% in the effective dielectric constant, which again is acceptable in 
the calculation of solvation energies. Finally, the site-potential representation 
{{ipa}) of [ini [H] exhibits the poorest convergence, particularly in the strong 
electric field case. Although the {ipa} entropy will eventually converge to the 
same value as that of the {puj} representation, the approximate {/i, represen- 
tation yields a more accurate free energy at practical iteration counts. 

Next, we turn to the convergence of the free energies with respect to the 
discretization of the orientation integrals. Figure [2] shows the relative error in 
the free energy for each orientation quadrature in Table A.l compared to the 
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3 4 5 6 7 8 9 10 5 10 100 105 110 

^ Distance [bohrs] 

(a) 4 A hard sphere (b) Capacitor with Eq = I V/A 

Figure 3: Convergence of density profiles with orientation quadrature for the two systems 
considered in Figure [l] The plotted site-densities are scaled by their corresponding bulk 
values, so that the profiles equilibrate at 1 far from the sphere / plates. Note that the 
densities at the lowest and highest quadratures are indistinguishable for the hard sphere, 
whereas the densities become similar to the fully-converged ones only around jmax ~ 10 for 
the strong-field capacitor. 



Euler(12) quadrature (taken to be the converged value) for the two systems 
considered above. The quadratures are sorted by jmax, the maximum degree 
of Wigner functions i^^^m^ which they are exact. (See Appendix A for 
details.) The relative error in the free energy decreases rapidly with quadrature 
size and plateaus ~ 10"^ at jmax ^ 7 for the hard sphere, limited by other 
discretization errors. For the highly polarized capacitor, higher quadratures 
are needed for the same level of accuracy, and the plateau occurs ~ 10~^ at 
jmax ^ 10. A reasonable choice for jmax for a typical system should therefore 
range from 7 to 10 depending on the strength of electric fields involved. 

Figure [3] shows the density profiles next to the hard sphere and the walls of 
the capacitor for various orientation quadratures. In the case of the hard sphere, 
the density profiles are virtually identical for all considered quadratures, as is 
expected given that the relative error in the free energy is ~ 10~^ even for 
the Octahedron group, one of the lowest quadratures considered with jmax = 
3. However, there are qualitative differences in the density profiles near the 
capacitor walls for jmax = 3 from the converged ones at jmax = 23 (Euler(12) 
quadrature), and the differences begin to disappear only around jmax — 10. At 
these field strengths, the orientation distribution is highly polarized (close to 
saturation) and hence requires a dense orientation quadrature to resolve. (The 
orientation distribution approaches a (5-function in the limit of infinite electric 
field.) 

4-.3. Accuracy of water Junctionals 

Finally, we turn to a comparison of the excess functionals for water suit- 
able for ah initio solvation methods. In particular, we focus on the scalar-EOS 
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Figure 4: Partial radial distributions (site-site correlation functions) for the scalar-EOS wa- 
ter functional (solid red lines), bonded- voids functional 27 (long-dashed green lines) and 
fittcd-correlations functional 1131 (short-dashed blue lines), compared to experimental pair 
correlations of water from Soper et al. |34j (black dotted lines). The position and location 
of the first goo peak for scalar-EOS and bonded-voids are in reasonable agreement with 
experiment, but the remaining structure resembles that of a close-packed hard sphere fluid 
rather than a tetrahedrally bonded one. The fitted-correlations functional is defined only at 
298 K and captures the features of the correlation functions by construction, but suffers from 
short-ranged artifacts due to the bandwidth-limited fitting procedure of | 13| . 



functional of section |3.2[ the bonded- voids functional [27] and the functional 
of Lisehner et al. |13j. The last functional is based on experimental corre- 
lations functions, which we will refer to as the 'fitted-correlations' functional. 
We perform all calculations in one-dimensional planar or radial grids, using the 
FluidlD sub-project of JDFTx [SD]. We use the Euler(20) orientation quadra- 
ture, with rict = 1 to exploit rotational symmetry in the transverse directions. 



(See Appendix B 



First we examine the pair correlation functions g^^ of the bulk fluid com- 
puted using the Ornstein-Zernike relation for the rigid-molecular fluid which 
may be written as 

h=(\- icNY^ici, (37) 

which is a matrix equation in Fourier space for each wave vector k. Here, 
hapik) is the Fourier transform of gap{r) — 1, Iap{k) = jo^kRa/s) is the intra- 
molecular structure factor with i?^^ being the distance between sites a and 
P within the molecule, N is the bulk number density of fluid molecules, and 
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Caii{k) is the Fourier transform of the direct correlation function Capir — r') = 
{-l/T)5'^F^^/5No,{r)5Np{r') evaluated in the limit of the uniform fluidQ 

The direct correlation functions are computed analytically in Fourier space 
for a set of wave vectors corresponding to the spherical Bessel function basis 



with 1024 basis functions and a radial extent r^ax = 64 bohrs (see Appendix 
[b| ), and the pair correlation functions are computed via (37) using numerical 
spherical Bessel transforms. Figure [4] compares the pair correlations for all three 
functionals under consideration compared against those obtained by Soper et al 
[34j from neutron diffraction data by empirical-potential structure refinement 
(EPSR). 

The scalar-EOS functional correctly captures the location and height of the 
first peak in gooif), but produces a secondary structure reminiscent of the 
close-packed coordination of the hard sphere fluid rather than the tetrahedral 
coordination exhibited by water. The split hydrogen peaks in the experimental 
data are fused into a single broader one with the same particle content. These are 
qualitatively the same features as the bonded-voids functional, but with slightly 
better agreement for the scalar-EOS functional. After all, the motivation for 
the scalar-EOS functional was to simplify the bonded-voids functional because 
it captured free energies of cavity formation reasonably despite not exhibiting 
features of tetrahedral correlation. The fitted-correlations functional reproduces 
some of the features of the experimental correlation functions by construction, 
but exhibits artifacts at short distances due to the bandwidth limitation in the 
fitting procedure for the correlations (and partly because it does not employ a 
hard sphere reference). 

Next, we examine the free energies of planar liquid- vapor interface for each 
functional. The calculations are performed on a one-dimensional planar grid of 
length 96 bohrs with 768 sample points and basis functions. For each temper- 
ature, the pressure is adjusted to the boiling point, which corresponds to equal 
chemical potentials and bulk grand free energy densities for the two phases. The 
hard sphere radius i?HS = 1.36 A for the scalar-EOS functional was determined 
by matching the interface energy obtained from such a calculation at 298 K to 
the experimental value for the surface tension 72.0 mN/mj^ Figure [s] compares 
the temperature dependence of this interface energy against experimental val- 
ues for the surface tension. The scalar-EOS functional captures the trend in the 
experimental data slightly better than the bonded-voids functional. 

The planar interface energies provide a means to calibrate the liquid func- 
tionals against experimental measurements, and the excellent agreement for the 
temperature dependence after adjusting the surface tension at one temperature 
is promising. However, the applicability of a functional for molecular solvation 



* The relation j37| may be generalized to mixtures of rigid-molecular fluids by replacing 
N with a diagonal matrix with the bulk number density of each component in the mixture, 
and setting /^^ = when a and /3 belong to different components of the mixture. 

^ The attraction range parameter au in the bonded-voids model 1271 and the smoothing 
parameter ro of the fitted-correlations model |13| were also fit to reproduce the surface tension 
at 298 K using similar calculations. 
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Figure 5: Energy of the planar vapor-liquid interface for the scalar-EOS and bonded-voids 
water functionals as a function of temperature, compared to the experimental values for surface 
tension 1351 . Both functionals fit the range parameter of the Lennard- Jones pair potential to 
the experimental surface tension at 298 K, and the scalar-EOS functional reproduces the 
temperature dependence more accurately than the bonded-voids one. (The fitted-correlations 
functional is omitted from this plot, since it is defined only at 298 K.) 
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Figure 6: Variation of the solvation energy of hard spheres that exclude the oxygen sites of 
water from their interior, with the radius of such spheres, compared to the SPC/E molecular 
dynamics results of |36| . The SPC/E model underestimates the bulk surface tension of water 
by 10% 1371 . and we have included a scaled version of the SPC/E data as a reasonable guess 
for real water. The scalar-EOS functional agrees with the bulk-scaled SPC/E data accurately, 
while the fitted-correlations functional systematically underestimates and the bonded-voids 
functional overestimates the free-energy of cavity formation. 
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Figure 7: Radial density profiles around hard spheres that exclude the oxygen sites of water 
from their interior, for spheres of radius 2, 4, 6, 8 and 10 A, compared to the SPC/E molecular 
dynamics results of [36| . The fitted-correlations functional misses the secondary peaks in the 
profiles, while the scalar-EOS and bonded-voids functional overestimate the contact density 
and the secondary structure, with best agreement provided by the scalar-EOS model. 

calculations depends on its ability to accurately describe the free energies re- 
quired to form cavities of molecular dimensions. A standard test case is the 
solvation free energy for microscopic hard spheres in the fluid. We compute the 
cavitation energies for hard spheres of radii R ranging from to 9 A, with exter- 
nal potentials Vo{r) = (1 Eh)0{R — r) and Vff (r) = that exclude the oxygen 
site of water from the interior of the spheres. The calculations are performed 
on a one-dimensional radial grid of extent r^ax = 64 bohrs (~ 34 A) with 512 
sample points and basis functions. 

Figure [6] compares the variation of the hard sphere solvation energy per sur- 
face area with sphere radius for all three functional with SPC/E molecular dy- 
namics estimates for the same from |36j. For large spheres, the surface curvature 
effects become negligible and the surface energy approaches the planar surface 
tension; whereas for small enough spheres the cavitation energy is proportional 
to the volume (A* = A^buikT x [A-KR^/i), so that A$/(47ri?2) cx R). Note that 
all three functionals agree perfectly with the molecular dynamics results in the 
small radius limit, and they all approach the bulk experimental surface tension 
in the large radius limit (after overshooting the bulk value in the bonded-voids 
case). However the SPC/E model underestimates the bulk surface tension to be 
65 mN/m f5T compared to the experimental value of 72 mN/m, and therefore 
the molecular dynamics results for the sphere solvation energies are also under- 
estimated by a similar amount for the larger spheres. Consequently, we include 
the molecular dynamics results scaled up by the ratio of experimental to SPC/E 
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Figure 8: Nonlinear dielectric response of the water functionals compared to SPC/E molec- 
ular dynamics results from | 38| . All three functionals provide essentially the same dielectric 
response, as this is determined by a competition between the molecular ideal gas entropy 
and the scaled mean field electrostatics. The minor differences arise from differences in the 
equations of state due to electrostriction (change of bulk-density in strong fields). The fitted- 
correlations functional has an unphysical instability accompanied by a rapid increase in density 
and drop in dielectric constant at an external field 2 V/A, due to an underestimation of 
the compressibility at high pressures by its polynomial excess free energy density model. 



surface tensions as a reasonable guess for the hard sphere cavitation energy of 
real water in Figure [6] (in addition to the unsealed values) The scalar-EOS 
functional significantly outperforms bonded-voids and fitted-correlations in its 
agreement with the bulk-scaled molecular dynamics results, and is the best can- 
didate for an accurate density functional description of cavitation energies in 
liquid water. 

We next examine the distribution of water around these spherical cavities 
of selected sizes in Figure [7] As expected from the results for the free energies, 
the density profiles of the scalar-EOS functional are in closest agreement with 
the SPC/E molecular dynamics results of [35]. The bonded- voids functional 
overestimates the structure in the liquid, which is expected since it also overes- 
timated the structure in the pair correlations (Figure |4]). The fitted-correlations 
functional severely underestimates the secondary structure in the density pro- 
files despite better qualitative agreement with the experimental pair correlation 
functions. 

Finally, we turn to the last key ingredient for a successful theory of solvation: 



The TIP4P/2005 pair potential for water captures the bulk surface tension much more 
accurately than SPC/E |37l . and it would be interesting to compare our density functional 
results to simulations of microscopic hard sphere solvation with that model. However, such 
results for TIP4P/2005 (analogous to |36| for SPC/E) have not yet been published to our 
knowledge. 
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nonlinear dielectric response. The typical electric fields in the vicinity of polar 
molecules are ~ V/A i.e. 10^° V/m, which corresponds to strong non-linearities 
and significant dielectric saturation. Here, we examine the nonlinear dielectric 
constant defined by e(£'o) = Eq/{Eq — inP), where Eq is a macroscopic exter- 
nally applied field and P is the corresponding bulk polarization density in the 
liquid. 

At equilibrium, a liquid in a macroscopic parallel-plate capacitor adopts uni- 
form density and polarization except for microscopic regions around the plates. 
The free energy of that capacitor is dominated by the bulk; the regions next to 
the plates only contribute via long-range interactions of the bound sheet-charge 
densities ±P in the liquid. Accounting for the interaction of these sheet charges 
with the external field, and with each other via the scaled mean-field Coulomb 
interaction, we can show that the effective free energy density minimized by the 
macroscopic capacitor at equilibrium is 

-E,-P + A,{T) — . 

Here, /ox is the excess-free energy density of the uniform fluid (which is de- 
termined entirely by the equation of state) and P = f ^^p^^ o P^oi is the 
polarization density, with Pmoi being the dipole moment of the fluid molecule 
in its reference orientation. We therefore minimize this free energy density on a 
planar grid with a single grid point to obtain the equilibrium P for each applied 
Eq, thereby avoiding the need for setting up a capacitor in a large simulation 
ceh. 

All three functionals considered here employ the same scaled mean-field elec- 
trostatic interaction constrained to produce the bulk dielectric response as pro- 
posed by Lischner et al. [IH]. The physics of dielectric saturation is captured 
by an interplay of this term with the entropy of the ideal gas of rigid molecules, 
which again is common to all three functionals. Consequently, their nonlin- 
ear dielectric response is very similar and compares quite well with the SPC/E 
molecular dynamics results |38j as shown in Figure U The minor differences 
between the functionals are due to the different uniform fiuid excess free energy 
densities (/ex) which correspond to different approximations to the equation of 
state of the fluid. The fitted-corrclations functional employs a polynomial model 
for /ox obtained from the bulk modulus and its pressure derivative at ambient 
conditions |T3j, which underestimates the bulk modulus at high compression. 
This causes the instability at high fields associated with a rapid increase in 
density, seen as a drop in the dielectric response at w 2 V/A in Figure [s] 

5. Conclusions 

We construct a general framework for the classical density-functional the- 
ory of rigid-molecular fluids that avoids the inversion problem associated with 
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site-density constraints by switching to the orientation density as the key vari- 
able. We show that the independent variables in previous solutions, such as 
ideal-gas effective site potentials, are compressed maximum-entropy representa- 
tions of the orientation density. We then motivate other representations with 
superior convergence properties which are variational approximations to the free 
energy. The self-representation, directly minimizing over the orientation den- 
sity {puj}, exhibits the fastest convergence for conjugate gradients minimization, 
but requires memory in proportion to the size of the quadrature for orientation 
integrals. The site-potential representation {ipa}, although exact in principle, is 
impractical due to poor convergence, particularly in the presence of strong elec- 
tric fields. We introduce the multipole representation {fj,e\ which exhibits com- 
parable convergence to the self-representation without the memory overhead, is 
effectively more accurate than {ipa} at practical iteration counts despite being 
a (variational) approximation, and enables efficient large-scale ab initio solva- 
tion in polar molecular fiuids within the framework of joint density-functional 
theory. 

We extend the algebraic formulation of electronic density-functional theory, 
DFT-I--I- [12], and present the discretization of our general framework and ex- 
cess functionals for practical calculations in a basis-independent manner. The 
methods developed in this paper form the basis for the fiuid sector of the open- 
source electronic density-functional theory software JDFTx [T3] , which provides 
a three-dimensional plane-wave basis implementation of this work. Addition- 



ally, a one dimensional version implementing the three basis sets of Appendix 
[B]is distributed as a sub-project of JDFTx j30], suitable for rapid prototyping 
and development of fluid functionals within this framework. 

In addition to the general framework for polar fluids, we construct a practical 
free energy functional for liquid water which improves on the accuracy of earlier 
functionals, the bonded-voids model |27| based on Wertheim perturbation and 
the fitted-correlations model |13j based on experimental correlation functions 
of water, while minimizing complexity and avoiding over-parametrization. We 
show that this 'scalar-EOS' functional accurately captures the key quantities of 
interest for ab initio solvation calculations: free energies for formation of mi- 
croscopic cavities in the fluid, and non-linear dielectric response. Within joint 
density-functional theory, the methods developed here provide an accurate and 
efficient description of solvent environments, thereby enabling a focused elec- 
tronic structure study of solvated biological and chemical systems of technolog- 
ical relevance. 

This work was supported as a part of the Energy Materials Center at Cornell 
(EMC^), an Energy Frontier Research Center funded by the U.S. Department of 
Energy, Office of Science, Office of Basic Energy Sciences under Award Number 
IdfcscDOOiOm 

Appendix A. Efficient quadratures for orientation integrals 

Efficient discretization of the orientation integrals is critical to the perfor- 
mance of any of the representations of Section |2.2| and determines the very 
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practicality of the Puj (self) representation. Here, we list efficient quadratures 
for discretizing integrals over w, J ^f{u}) Wif{u}i). 

The simplest approach is to label orientations by ZYZ-Euler angles to = 
(a, P, 7) and use the outer product of a Gauss-Legendre quadrature for /3 G [0, tt] 
and Gauss- Fourier quadratures for the periodic a, 7 € [0, 27r). More efficient 
quadratures may be constructed as an outer product using the S2 x Si structure 
of 50(3), or by working directly on S0{3) without an outer product structure 
[39]. 

In |39j . quadratures on 50(3) are optimized to minimize the RMS error in 
the integrals of all D^^j^^{uj) up to some jmax- We focus on quadratures that 
are exact up to some jmax, 



E 



W^.^^;,™.(^.) = E^»<i™2(/30e^^'"^"-+'"^"'^ = Sjo (A.l) 



for all \mi\, \m2\ < i < Jmax: and can be optimized further using the symmetry 
of the molecule at hand. For simplicity, we only consider Z„ symmetry about 
a single axis, chosen to be the 2;-axis of the molecule frame without loss of 
generality. The quadratures considered then fall into 3 classes: 

(a) Symmetry groups of Platonic solids [32] 

(b) Outer products of a spherical j-design [JD] on §2(0, /3) with a uniform 
quadrature on §1(7) 

(c) Outer product quadrature on all 3 Euler angles a, f3 and 7. 

Each of these these classes consists of uniformly spaced nodes of equal weights 
in 7 for each (a,/3). Grouping the nodes as (afc,/3fc,7fc + 2n7r/n^) for n € 



0, . . . , — 1 with total weight Wk for each group, (A.l I can be reduced to 

E Wkdi^.^My^"''""^"''''^ = (A.2) 

k 

for all |mi|, \m2\ < j < jmax such that TO2 is a multiple of n-y. Therefore if 
n-y > jmax (which is the case for all but the Icosahedron rotation group), (A.2 1 
further simplifies to 

^ WkY^, {/3k , ttfe) = S,Q (A.3) 

k 

for all |m| < j < jmax using the relations of D^^^^ to the spherical harmonics. 
A spherical jmax-design is a set of points on the unit sphere that satisfies 



(A.3) with uniform weights Wk, and hence it yields an 5*0(3) quadrature exact 
to jmax when combined with a uniform quadrature with jmax + 1 nodes on 
Si (7). We use the spherical designs with the smallest number of nodes for each 
7 < jmax < 21 tabulated in ^40 to form the quadratures of class (b). The 
quadratures of lower order reduce to class (a), specifically the rotation groups 
of the Tetrahedron at jmax = 2, Octahedron at jmax — 3 and Icosahedron at 

jmax — 5. 

The Gauss-Legendre quadrature with nodes on cos/3 G [^Ijl] is exact 
for the integration of all polynomials up to order 2np — 1. The outer product 
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of this with a uniform quadrature with 2np nodes on a £ [0, 2tt) satisfies (A. 3 1 



for jmax = 2n^ — 1, and hence also ( A.l I to that order when combined with 2n 



^/3 



uniform samples on 7. 

Finally the reduction by Z„ symmetry about the z-axis in the molecule 
frame amounts to replacing 81(7) with §i/Z„. This is achieved by a uniform 
sampling of [n-y/n] points on 7 g [0, 27r/n), which retains the exactness to jmax 
for functions with this symmetry with a reduction of up to n in the number of 
nodes required. 

The accuracy of these quadratures for the classical density functional theory 
of rigid molecules is explored in section |4.2| The quadratures considered there 
are listed in Table along with their jmax, the number of nodes for sampling 
SO{3) /Zn in general and S'0(3)/Z2 in particular, which is the case relevant 
for water. Note that the Euler quadrature with np — 3 needs almost twice as 
many nodes as the Icosahedron group for the same jmax = 5, but the relative 
inefficiency of the Euler quadratures decreases with j^ax and becomes less than 
1% between the np — 11 Euler quadrature and the 21-design at jmax = 21. 

Appendix B. One-dimensional discretization for special geometries 



The discretization of three-dimensional space according to Section [4~1] along 
with the orientation quadratures of [Appendix A provide a practical route to 



computations with the rigid-molecular classical density functional framework 
of Section [2] in arbitrary geometries and basis sets. However, the development 
and testing of new excess functionals for liquids primarily require calculations 
in high-symmetry configurations. Here, we detail the formulation of highly- 
efficient discretizations of planar, cylindrical and spherical geometries on a one- 
dimensional grid, which allow for the rapid prototyping of excess functionals 



employed in Section 4.3 and [27] . 

The discretization of space is performed in the framework of Section |4.1[ 
but with special basis sets exploiting the symmetry. The three geometries we 
consider here are 

1. Planar, where all spatial dependence is along z, 

2. Cylindrical, with dependence only on the distance from the z-axis p, and 

3. Spherical, with dependence only on distance from origin r. 

Each of these geometries require only a one-dimensional discretization. For the 
planar geometry, we impose mirror-symmetry boundary conditions at the ends 
of the grid, and pick a basis of cosines and a corresponding quadrature grid 
suited for the Discrete Cosine Transform [41] . For the spherical and cylindrical 
geometries, we impose Neumann boundary conditions at some maximum radius, 
and choose a finite basis of spherical and cylindrical Bessel functions respectively, 
along with a quadrature grid suited for the Discrete Bessel Transform [32] 



^^The Discrete Bessel Transform of '421 is based on Dirichlet boundary conditions; the 
extension of tliat approach to Neumann boundary conditions is straightforward, and the 
results are summarized in Table IB. 21 
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Jmax 


Number of quadrature nodes for 




50(3)/Z„ 


50(3)/Z2 


Tetrahedron 


2 


4 X [3/n] 


8 


Octahedron 


3 


6 X [4/n] 


12 


Icosahedron 


5 


12 X [5/n] 


36 


7-design 


7 


24 X [8/n] 


96 


8-design 


8 


36 X [9/n] 


180 


9-design 


9 


48 X [10/n] 


240 


10-design 


10 


60 X \ll/n\ 


360 


11-design 


11 


70 X [12/n] 


420 


12-design 


12 


84 X [13/n] 


588 


13- design 


13 


94 X [14/n] 


658 


14-design 


14 


108 X [15/n] 


864 


15-design 


15 


120 X [16/n] 


960 


16- design 


16 


144 X [17/n] 


1296 


17-design 


17 


156 X [18/n] 


1404 


18-design 


18 


180 X [19/n] 


1800 


19- design 


19 


204 X [20/n] 


2040 


20-design 


20 


216 X [21/n] 


2376 


21-design 


21 


240 X [22/n] 


2640 


Euler(n^) 


2np - 1 


2n| X \2np/n\ 


2n3 



Table A.l: List of explored quadratures, their degree of exactness Jmax, and the number of 
nodes in sampling SO{3)/Z,„. The Euler angles corresponding to the platonic solid rotation 
groups are listed in 1391 . The j'-designs are constructed as an outer product of the spherical j- 
designs with fewest points for each j from [40] used for (o, f}) with \{j + \ uniform samples 
on 7 £ [0, 2-K jn). Each Euler(?i^) quadrature is an outer product of a n^-point Gauss-Legendre 
quadrature on cos /3 £ [—1, 1], a uniform 2n^-point quadrature on a £ [0, 27r), and a uniform 
[2n^/n] -point quadrature on 7 G [0, 27r/n). 
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Planar 


Cylindrical 


Spherical 


Coordinate 
System 




{p, <l>,z) 




Symmetry 


f{r) -> .f{z) 


^ f{p) 


f{r) ^ fir) 


Boundary 
conditions 


f{0) = f{L)=0 


/'(Pmax) = 


/'(n„ax) = 


RpQiQ h • { r\ 


Wi cos(Giz), 
Gi = in/L, 

"^^ ~ (l+5io)i 


t2;jJo(Gip), 


Gi — yz/^inax; 


Quadrature 
grid {fj} 








T 


W,;COS {{j + 5)7r|) 






Jij 


WjCos{{j + Dtt^) , 




J ~ !/sif(^j + i) 




WiSi'i 




-G^WiSi'i 








Si'i J 47rr2drg(r)jo(G,r) 



Table B.2: Definition of the basis functions for the high-symmetry geometries - planar, cyhn- 
drical and spherical - with one-dimensional discretizations of sample count S, and matrix 
elements of the operators of Section |4.1| for each of these basis sets. The basis functions are 
labeled by i = 0, 1, ■ ■ • ,5 — 1 for each basis set, and Xi, Xi, Yi and j/i, are the i^^ roots of 
Jo(x), jo(x), Jq{x) and j'gix) respectively, with Yq = yo = 0. The quadrature grid has the 
same number of points S as the basis size, and are labeled by j = 0, 1, • • • , S — 1. 



The definition of the basis functions, quadrature grid and the matrix elements 



for the operators of Section 4.1 are summarized in Table B.2 



All three basis sets are derived from the eigenfunctions of the three-dimensional 
Laplace equation in various geometries, and are therefore intricately linked to 
the three-dimensional plane- wave basis: the basis functions are indexed by Gi, 
the magnitude of the corresponding plane-wave momentum. Consequently, the 
Laplacian and convolutions by spherical functions are diagonal in these basis 
sets as well, as indicated in Table B.2 The transform operators I and J reduce 
to the 'DCT type III' and 'DCT Type IT fast Fourier transforms [43] respec- 
tively in the planar geometry (or 'IDCT' and 'DCT' in the notation of [4T]): 
the cylindrical and spherical transforms lack an analogous ©(Slogs') method 
and are implemented as matrix-vector multiplies with a precomputed Bessel 
function matrix. 

The basis-independent discretization of the scalar-EOS excess functional 
(35), and site-density excess functionals in general, carries over to the planar, 
cylindrical and spherical geometries without modification. The discretization 
of the rigid-molecular ideal gas free energy and the generation of site-densities 
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from independent variables carries over unmodified for the planar geometry, but 
is slightly complicated for the cylindrical and spherical geometries by the fact 
that the translation operator breaks the symmetry of the basis set and does not 
have a one-dimensional representation. 



We can however compute the site-densities using (34) and the orientation 



density in the site-potential representation using (31 1 for these basis sets as 
well, with minor modifications to the translation operators in those equations. 
First, we pick a covariant reference orientation for the molecule, (relative to 
the local coordinate frame (p, (j), z) or (f, 9, (f))), so that Puj{r) is invariant under 
the cylindrical or spherical symmetry for each w and permits a one-dimensional 
representation!^ Consequently, the translations involved in (34 1 and (31 ) would 
be relative to the local coordinate frame as well, and hence position-dependent; 
we therefore need to generalize the translation operators 7s to 'warp' operators 
7s(f?) defined by 7s(f)/(r) = f[r + d{r)). It can be shown that the expressions 



of Section 4.1 remain valid without modification upon this generalization. 

The translation operator for the planar basis is a simple one-dimensional 
restriction of its three-dimensional counterpart, and it generalizes to 

rs(p)/(p) - / (^\/(p + 5-/5)2 + (a-<^)2j (B.l) 

for the cylindrical basis with f{p) = /(2/3,uax ^ p) for p > Pmaxi and 

Tsir)f{r) = / (v/r2 + a2 + 2ra.f) (B.2) 

for the spherical basis with f{r) = /(2rmax — ioi r > rmaxj^ We could 
compute the matrix elements of these operators in the Bessel basis and apply 
the translation as a dense-matrix multiply in basis space, but those suffer from 
Nyquist frequency ringing problems similar to their three-dimensional counter- 
parts. Instead, we compute these operators in real space using approximate 
sampling operators 5s(f) based on constant or linear-spline interpolation which 
preserve non-negativity of scalar fields. 

The results for the scalar-EOS water functional in Section 14.31 and the 
bonded-voids water functional in |27j were computed using the discretization 
scheme of Section [41] in the planar and spherical bases, with the warp operator 
S computed using linear-spline interpolation as discussed above. The planar 
and spherical bases have an additional rotational symmetry about the local z 
and f axes respectively at any point in space which renders the integral over 
Euler angle a trivial, so that a quadrature on 82(7,/?) with no a sampling suf- 
fices; the one-dimensional calculations employ this additional optimization by 



^^If wc used an invariajit reference orientation as in the three-dimensional case, Pui{r) would 
be covariant under the symmetry, so that the spatial dependence of p^j (r) for each ui would 
not be cylindrically or spherically symmetric, and would therefore lack a one-dimensional 
representation. 

The covariant reference frame ensures that a ■ p and a • cj> depend only on p (and not (f) 
and z), and that a - f depends only on r. 
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using the Euler(7i^) quadratures of Appendix A but with Ua — 1 irrespective 

of 71^. 
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